library("dplyr")
library("raster")
library("sf")
library("data.table")
library("datawizard")
library("zoo")
library("doParallel")
library("foreach")
library("parallel")
library("ncdf4")
library("SPEI")
library("ggplot2")
library("foreign")
library("tidyr")
library("terra")
library("stars")
library("yarrr")
library("wesanderson")
#library("ClimateOperators")
library("rlist")
library("writexl")
library("readxl")
library("stringr")
library("tidyr")
library("miceadds")
library("stargazer") #Display table on latex format
library("foreign")
library("dplyr") #data management
library("tidyverse") #data management
library("tibble") # tibbles
library("caret")
library("naniar")
library("mapproj")
library("ggmap")
library("RgoogleMaps")
library("mapview")
library("leaflet")
library("haven")
library("rnaturalearth")
library("reshape2")
library("sjlabelled")
library("RColorBrewer")Tutorial Functions average distance and average density
1 Tutorial to calculate average distance and point density statistics
This guide provides a step-by-step approach to compute the average distance and point density statistics for a target set of polygons/points coordinates in R. It uses own functions wrappers: ave_dist_2shp() and ave_point_density_2shp(). The inputs are shapefile files loaded as sf_dataframe R objects and the returning outputs are dataframe R objects.
2 Overview of Steps
We will go through the following steps:
- Load the data
- Prepare the data
- Compute the average distance to rivers for each
ID(village, adm-div, etc.) - Compute the average point density of buildings for each
ID(village, adm-div, etc.)
3 Loading packages
All packages loaded
4 Setting up working directory and path to wrapper functions
Working directory set up as default with folder’s location.
# Otherwise use
#setwd("C:/Users/X/OneDrive - Food and Agriculture Organization/PROJECT11_GEOINDICATORS_REPOSITORY")
rm(list = ls()) # remove existing objects# load wrapper functions
source(file.path(#"web_tutorial",
"functions.R"))
path_to_project <- file.path("..", "..",
"data", "data_project", "benin",
"PROJECT11_GEOINDICATORS_REPOSITORY")5 Calculating average distance to rivers for each polygon ID
6 Loading inputs
All input have to be set up as EPSG4326 or WGS84 as coordinate system
6.0.1 Shapefile extension for the target country (Benin)
This uses GADM database (see GADM)
path_gadm3 <- file.path(path_to_project,
"INPUT",
"Benin_shp","Benin-ADMlevels","Shapefile",
"BEN_shp","gadm41_BEN_3.shp")
gadm3 <- st_as_sf(st_read(path_gadm3)) |>
dplyr::select(-NL_NAME_1, -NL_NAME_2, -NL_NAME_3,
-VARNAME_3, -CC_3, -HASC_3,-TYPE_3,
-ENGTYPE_3,-COUNTRY)gadm3 <- st_transform(gadm3, 4326) #defining as EPSG4326 or WGS84
gadm3Simple feature collection with 546 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 0.774345 ymin: 6.23491 xmax: 3.851701 ymax: 12.41835
Geodetic CRS: WGS 84
First 10 features:
GID_3 GID_0 GID_1 NAME_1 GID_2 NAME_2 NAME_3
1 BEN.1.1.1_1 BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Banikoara
2 BEN.1.1.2_1 BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Founougo
3 BEN.1.1.3_1 BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Gomparou
4 BEN.1.1.4_1 BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Goumori
5 BEN.1.1.5_1 BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Kokey
6 BEN.1.1.6_1 BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Kokiborou*
7 BEN.1.1.7_1 BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Ounet
8 BEN.1.1.8_1 BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Somp?r?kou
9 BEN.1.1.9_1 BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Soroko
10 BEN.1.1.10_1 BEN BEN.1_1 Alibori BEN.1.1_1 Banikoara Toura
geometry
1 MULTIPOLYGON (((2.415323 11...
2 MULTIPOLYGON (((2.346528 11...
3 MULTIPOLYGON (((2.464541 11...
4 MULTIPOLYGON (((2.070032 11...
5 MULTIPOLYGON (((2.802918 11...
6 MULTIPOLYGON (((2.26218 11....
7 MULTIPOLYGON (((2.405739 11...
8 MULTIPOLYGON (((2.463629 11...
9 MULTIPOLYGON (((2.338276 11...
10 MULTIPOLYGON (((2.317923 11...
6.0.2 Shapefile target country (Benin)
This file is the shapefile with the polygon’s structure. Each ID is enumerated in column id_gns. The original file with X, Y villages’ coordinates can be found here Villages Benin. The file was modified and non-intersecting buffers (using Voronoi’s structure) of 2km were constructed using the villages’ coordinates.
path_polygons_gns <- file.path(path_to_project,
"INPUT","QGIS_map_voronoi",
"villages_voronoibuffer_adm4_gns.shp")
polygons_gns <- st_as_sf(st_read(path_polygons_gns))polygons_gns <- st_transform(polygons_gns, 4326) #defining as EPSG4326 or WGS84
polygons_gnsSimple feature collection with 3599 features and 18 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 0.7645527 ymin: 6.222888 xmax: 3.839671 ymax: 12.42317
Geodetic CRS: WGS 84
First 10 features:
ID gid_0 gid_1 gid_2 gid_3 name_1 name_2 name_3
1 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
2 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
3 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
4 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
5 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
6 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
7 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
8 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
9 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
10 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
type_3 engtype_3 id_gns full_name full_nm_nd sort_name desig_cd
1 Arrondissement Borough 2058 Gano Gano GANO PPL
2 Arrondissement Borough 1426 Kokabo Kokabo KOKABO PPL
3 Arrondissement Borough 973 Wagui Wagui WAGUI PPL
4 Arrondissement Borough 3392 Atabénou Atabenou ATABENOU PPL
5 Arrondissement Borough 2726 Komon Komon KOMON PPL
6 Arrondissement Borough 2263 Dérou Derou DEROU PPL
7 Arrondissement Borough 3386 Banikoara Banikoara BANIKOARA PPL
8 Arrondissement Borough 453 Toké-Banta Toke-Banta TOKEBANTA PPL
9 Arrondissement Borough 1873 Gomparou Gomparou GOMPAROU PPL
10 Arrondissement Borough 581 Sonworé Sonwore SONWORE PPL
adm1 lat_y lon_x geometry
1 BJ-000 11.24685 2.414291 MULTIPOLYGON (((2.423529 11...
2 BJ-000 11.28007 2.427724 MULTIPOLYGON (((2.434365 11...
3 BJ-000 11.28146 2.400951 MULTIPOLYGON (((2.392831 11...
4 BJ-000 11.28761 2.381083 MULTIPOLYGON (((2.392831 11...
5 BJ-000 11.28951 2.423543 MULTIPOLYGON (((2.426221 11...
6 BJ-000 11.29832 2.401679 MULTIPOLYGON (((2.413442 11...
7 BJ-AL 11.29845 2.438561 MULTIPOLYGON (((2.447332 11...
8 BJ-000 11.30770 2.414829 MULTIPOLYGON (((2.430931 11...
9 BJ-AL 11.33134 2.441141 MULTIPOLYGON (((2.437974 11...
10 BJ-000 11.35827 2.422149 MULTIPOLYGON (((2.440276 11...
6.0.3 Shapefile shape_input
We will use the shapefile of African Rivers from FAO as our shape_input (see FAO rivers1 and here FAO rivers2)
path_rivers_fao <- file.path(path_to_project, "INPUT","FAO_GIS","rivers_africa","rivers_africa_37333.shp")
rivers_fao <- st_as_sf(st_read(path_rivers_fao))rivers_fao <- st_transform(rivers_fao, 4326) #defining as EPSG4326 or WGS84
rivers_faoSimple feature collection with 185730 features and 18 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: -17.23958 ymin: -34.76458 xmax: 54.35 ymax: 37.21875
Geodetic CRS: WGS 84
First 10 features:
FID_af_str ARCID FROM_NODE TO_NODE FID_sub_ba SUB_BAS MAJ_BAS
1 0 1 2 1 32890 201591 7020
2 1 2 5 4 32890 201591 7020
3 2 3 6 7 32890 201591 7020
4 3 4 7 3 32890 201591 7020
5 4 5 8 6 32890 201591 7020
6 5 6 10 8 32890 201591 7020
7 6 7 11 9 32890 201591 7020
8 7 8 13 8 32890 201591 7020
9 8 9 14 6 32890 201591 7020
10 9 10 15 7 32890 201591 7020
MAJ_NAME SUB_NAME MAJ_AREA LEGEND SUBBAS_ID
1 Mediterranean South Coast Algerian east coast 558292 20 7201591
2 Mediterranean South Coast Algerian east coast 558292 20 7201591
3 Mediterranean South Coast Algerian east coast 558292 20 7201591
4 Mediterranean South Coast Algerian east coast 558292 20 7201591
5 Mediterranean South Coast Algerian east coast 558292 20 7201591
6 Mediterranean South Coast Algerian east coast 558292 20 7201591
7 Mediterranean South Coast Algerian east coast 558292 20 7201591
8 Mediterranean South Coast Algerian east coast 558292 20 7201591
9 Mediterranean South Coast Algerian east coast 558292 20 7201591
10 Mediterranean South Coast Algerian east coast 558292 20 7201591
TOBAS_ID Strahler A_Strahler RASTERVA_2 RASTERVA_1 Regime
1 -999 1 8 609 5202 P
2 -999 1 8 460 2537 P
3 -999 2 7 518 97325 P
4 -999 3 6 478 218144 P
5 -999 2 7 518 41565 P
6 -999 1 8 518 8568 P
7 -999 1 8 413 9945 P
8 -999 1 8 638 5032 P
9 -999 1 8 648 6162 P
10 -999 2 7 478 106369 P
geometry
1 MULTILINESTRING ((9.220833 ...
2 MULTILINESTRING ((9.945833 ...
3 MULTILINESTRING ((9.65625 3...
4 MULTILINESTRING ((9.735417 ...
5 MULTILINESTRING ((9.647917 ...
6 MULTILINESTRING ((9.614583 ...
7 MULTILINESTRING ((10.075 37...
8 MULTILINESTRING ((9.479167 ...
9 MULTILINESTRING ((9.2375 37...
10 MULTILINESTRING ((9.68125 3...
As the shapefile of rivers corresponds to all Africa, we will crop it to country’s area using the extension from the object gamd, which has the extension of Benin. It will be enlarged a bit in case to guarantee it does contain all country’s borders. Notice also that for using the function raster::crop, we need to define the object as an Spatial one, do the cropping and define it again as sf_dataframe, for which, we need to define it again with the same coordinate system WGS84 (otherwise our function ave_point_density_2shp will trew an error).
ext_benin<-raster::extent(gadm3)
ext_benin@xmin<-ext_benin@xmin-0.5
ext_benin@xmax<-ext_benin@xmax+0.5
ext_benin@ymin<-ext_benin@ymin-0.5
ext_benin@ymax<-ext_benin@ymax+0.5
#cropping
rivers_fao <- st_as_sf(raster::crop(as(rivers_fao, "Spatial"),ext_benin))
rivers_fao <- st_transform(rivers_fao, 4326) #defining as EPSG4326 or WGS84 again as I only defined as sf dataframe before
rivers_faoSimple feature collection with 1796 features and 18 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 0.274345 ymin: 5.78125 xmax: 4.351701 ymax: 12.91835
Geodetic CRS: WGS 84
First 10 features:
FID_af_str ARCID FROM_NODE TO_NODE FID_sub_ba SUB_BAS MAJ_BAS MAJ_NAME
1 85827 85828 89460 89298 32706 20872 7002 Niger
2 85828 85829 89461 89407 32705 20871 7002 Niger
3 85829 85830 89348 89462 32705 20871 7002 Niger
4 85871 85872 89334 89502 32706 20872 7002 Niger
5 85872 85873 89502 89460 32706 20872 7002 Niger
6 85918 85919 89264 89462 32705 20871 7002 Niger
7 85930 85931 89569 89381 32706 20872 7002 Niger
8 85931 85932 89070 89571 32692 20843 7002 Niger
9 85979 85980 88722 89617 32701 20854 7002 Niger
10 85980 85981 89553 89617 32701 20854 7002 Niger
SUB_NAME MAJ_AREA LEGEND SUBBAS_ID TOBAS_ID Strahler A_Strahler
1 Faga 2136941 2 7020872 7020871 4 5
2 Niger 9 2136941 2 7020871 7020859 1 8
3 Niger 9 2136941 2 7020871 7020859 6 3
4 Faga 2136941 2 7020872 7020871 1 8
5 Faga 2136941 2 7020872 7020871 4 5
6 Niger 9 2136941 2 7020871 7020859 4 5
7 Faga 2136941 2 7020872 7020871 1 8
8 Sokoto 1 2136941 2 7020843 7020841 1 8
9 Dallol Maouri 2136941 2 7020854 7020853 2 7
10 Dallol Maouri 2136941 2 7020854 7020853 1 8
RASTERVA_2 RASTERVA_1 Regime geometry
1 182 907940 P MULTILINESTRING ((0.5229167...
2 421 10060 I MULTILINESTRING ((2.102083 ...
3 55 66217700 P MULTILINESTRING ((2.376455 ...
4 667 5300 I MULTILINESTRING ((0.4721039...
5 674 894820 P MULTILINESTRING ((0.4979167...
6 158 1462580 P MULTILINESTRING ((2.256649 ...
7 262 5870 I MULTILINESTRING ((0.76875 1...
8 280 6090 I MULTILINESTRING ((4.34375 1...
9 705 178340 I MULTILINESTRING ((3.49375 1...
10 586 5100 I MULTILINESTRING ((3.514583 ...
7 Calculating average distance to rivers for each ID polygon in target
This wrapper function calculates the average distance (in m) of each ID element in target to the object shape_input. It uses distanceFromPoints which works better, faster than st_distance function. For the moment, the function works better for UTM coordinate system.
country_ext defines the whole extension of the country, including country borders. Here, we use gadm3 is our country’s extension for Benin.
target is the shapefile (adm-div such as communes, villages or HH coordinates) for which the average distance will be calculated. We use polygons_gns as our target polygon for which we will calculate the distance.
shape_input is the shapefile for which a raster of distances will be calculated (eg. roads, rivers, etc.). In our case, we use rivers_fao as the object from which the distances will be calculated.
UTM corresponds to the country’s location in the UTM coordinate system. As for calculating distances it is more convenient to use meters instead of degrees, and the coordinate system needs to be transformed into the UTM of the country, which is more precise around the equator. We define UTM=31, which corresponds to Benin. Otherwise, it uses default at the global level (3857) (under development).
res_raster is the resolution in meters for the raster of distances that will be created. A higher resolution increase the time of processing. We define a resolution of 10000 which corresponds to 10 km. Internally, it will create a raster of distances of 10 km, and each pixel size of 10kmX10km.
name_dist name of the average distance variable constructed. It will be created as ave_dist_NAME with name_dist=“NAME” giving the name of the distance variable that will be created. In our case, ave_dist_rivers.
ave_dist_rivers <- ave_dist_2shp(country_ext=gadm3,target=polygons_gns,shape_input=rivers_fao,UTM=31,res_raster=10000,name_dist="rivers")[1] "UTM is: 31"
ave_dist_rivers[1:10, ] #print first 10 elements ID gid_0 gid_1 gid_2 gid_3 name_1 name_2 name_3
1 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
2 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
3 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
4 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
5 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
6 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
7 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
8 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
9 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
10 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
type_3 engtype_3 id_gns full_name full_nm_nd sort_name desig_cd
1 Arrondissement Borough 2058 Gano Gano GANO PPL
2 Arrondissement Borough 1426 Kokabo Kokabo KOKABO PPL
3 Arrondissement Borough 973 Wagui Wagui WAGUI PPL
4 Arrondissement Borough 3392 Atabénou Atabenou ATABENOU PPL
5 Arrondissement Borough 2726 Komon Komon KOMON PPL
6 Arrondissement Borough 2263 Dérou Derou DEROU PPL
7 Arrondissement Borough 3386 Banikoara Banikoara BANIKOARA PPL
8 Arrondissement Borough 453 Toké-Banta Toke-Banta TOKEBANTA PPL
9 Arrondissement Borough 1873 Gomparou Gomparou GOMPAROU PPL
10 Arrondissement Borough 581 Sonworé Sonwore SONWORE PPL
adm1 lat_y lon_x ave_dist_rivers
1 BJ-000 11.24685 2.414291 4885.229
2 BJ-000 11.28007 2.427724 4885.229
3 BJ-000 11.28146 2.400951 4885.229
4 BJ-000 11.28761 2.381083 4885.229
5 BJ-000 11.28951 2.423543 4885.229
6 BJ-000 11.29832 2.401679 4885.229
7 BJ-AL 11.29845 2.438561 5402.293
8 BJ-000 11.30770 2.414829 4885.229
9 BJ-AL 11.33134 2.441141 6628.781
10 BJ-000 11.35827 2.422149 6761.864
8 Calculating average density to point-buildings for each ID polygon in target
This wrapper function calculates the average density of points for each ID element of targetto the object shape_input .
9 Using Buildings dataset as shape_input
While some of the objects are already defined as in the previous function (gadm3 and polygons_gns), we will use as shape_input a shapefile that contains the points with buildings located in the country (see OCHA Datahum). The file comes in zip format as it is opened, transformed as an sf_dataframe
path2buildings <- file.path(path_to_project,
"INPUT", "Benin_OCRA","Buildings","hotosm_ben_buildings_polygons_shp.zip") #using zip
# Extract the ZIP file and Find the file among the extracted files. There is only one
output_directory_buildings <- tempdir() # Using a temporary directory
unzip(path2buildings, exdir = output_directory_buildings)
extracted_files <- list.files(output_directory_buildings, full.names = TRUE)
buildings_file_path <- extracted_files[grep("\\.shp$", extracted_files)]
#reading file
buildings_ocra <- st_as_sf(st_read(buildings_file_path))buildings_ocra <- st_transform(buildings_ocra, 4326) #defining as EPSG4326 or WGS84
buildings_ocraSimple feature collection with 813093 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 0.7848944 ymin: 6.236586 xmax: 3.836711 ymax: 12.40539
Geodetic CRS: WGS 84
First 10 features:
osm_id addrcity office name buildingle addrhousen
1 81452793 <NA> <NA> Tri Postal <NA> <NA>
2 81766299 <NA> <NA> <NA> 8 <NA>
3 116969223 <NA> <NA> Site Pompe <NA> <NA>
4 116969224 <NA> <NA> Station SONAB <NA> <NA>
5 116969225 <NA> <NA> Site Pompe <NA> <NA>
6 116985473 <NA> <NA> SBEE <NA> <NA>
7 116985474 <NA> <NA> Hôtel de Ville de Bassila <NA> <NA>
8 116985476 <NA> <NA> Hôtel de Ville de Bassila <NA> <NA>
9 167008127 <NA> <NA> <NA> <NA> <NA>
10 167008128 <NA> <NA> <NA> <NA> <NA>
addrfull building addrstreet source buildingma
1 <NA> yes <NA> <NA> <NA>
2 <NA> yes <NA> <NA> <NA>
3 <NA> yes <NA> <NA> <NA>
4 <NA> yes <NA> <NA> <NA>
5 <NA> yes <NA> <NA> <NA>
6 <NA> yes <NA> <NA> <NA>
7 <NA> yes <NA> <NA> <NA>
8 <NA> yes <NA> <NA> <NA>
9 <NA> yes <NA> <NA> <NA>
10 <NA> yes <NA> <NA> <NA>
geometry
1 MULTIPOLYGON (((2.386848 6....
2 MULTIPOLYGON (((2.386307 6....
3 MULTIPOLYGON (((1.666625 9....
4 MULTIPOLYGON (((1.667526 9....
5 MULTIPOLYGON (((1.6682 9.00...
6 MULTIPOLYGON (((1.667723 8....
7 MULTIPOLYGON (((1.667386 8....
8 MULTIPOLYGON (((1.667395 8....
9 MULTIPOLYGON (((2.340467 6....
10 MULTIPOLYGON (((2.340359 6....
Below we describe the inputs for the function:
country_ext defines the whole extension of the country, including country borders. Here, we use gadm3 is our country’s extension for Benin.
target is the shapefile (adm-div such as communes, villages or HH coordinates) for which the average point density will be calculated. We use polygons_gns as our target polygon for which we will calculate the point density.
shape_input is the shapefile from which to calculate points density (eg. buildings, roads, rivers, etc.) In our case, it will be our building’s object buildings_ocra.
res_raster is the resolution in degrees for the raster of distances that will be created and applied to `shape_input’ (0.1 could be equivalent to 10km or 0.1*100000=10000m around the Equator). A higher resolution increases the time of processing.
name_dist name of the average density variable constructed. It will be created as ave_den_NAME with name_density=“NAME In our case, ave_den_allbuildings.
ave_den_allbuildings <-ave_point_density_2shp(country_ext=gadm3,
target=polygons_gns,
shape_input=buildings_ocra,
res_raster=0.1,
name_density="allbuildings")
ave_den_allbuildings[1:10, ] #print first 10 elements ID gid_0 gid_1 gid_2 gid_3 name_1 name_2 name_3
1 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
2 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
3 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
4 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
5 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
6 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
7 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
8 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
9 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
10 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
type_3 engtype_3 id_gns full_name full_nm_nd sort_name desig_cd
1 Arrondissement Borough 2058 Gano Gano GANO PPL
2 Arrondissement Borough 1426 Kokabo Kokabo KOKABO PPL
3 Arrondissement Borough 973 Wagui Wagui WAGUI PPL
4 Arrondissement Borough 3392 Atabénou Atabenou ATABENOU PPL
5 Arrondissement Borough 2726 Komon Komon KOMON PPL
6 Arrondissement Borough 2263 Dérou Derou DEROU PPL
7 Arrondissement Borough 3386 Banikoara Banikoara BANIKOARA PPL
8 Arrondissement Borough 453 Toké-Banta Toke-Banta TOKEBANTA PPL
9 Arrondissement Borough 1873 Gomparou Gomparou GOMPAROU PPL
10 Arrondissement Borough 581 Sonworé Sonwore SONWORE PPL
adm1 lat_y lon_x ave_den_allbuildings
1 BJ-000 11.24685 2.414291 804.0000
2 BJ-000 11.28007 2.427724 804.0000
3 BJ-000 11.28146 2.400951 804.0000
4 BJ-000 11.28761 2.381083 654.0000
5 BJ-000 11.28951 2.423543 804.0000
6 BJ-000 11.29832 2.401679 804.0000
7 BJ-AL 11.29845 2.438561 804.0000
8 BJ-000 11.30770 2.414829 595.8421
9 BJ-AL 11.33134 2.441141 92.1000
10 BJ-000 11.35827 2.422149 13.0000
10 Using points rather than non-intersecting buffers
Testing function for points:
path_points_gns <- file.path(path_to_project,
"INPUT","QGIS_map_voronoi","villages_points_adm4_gns.shp")
points_gns <- st_as_sf(st_read(path_points_gns))points_gns <- st_transform(points_gns, 4326) #defining as EPSG4326 or WGS84
points_gnsSimple feature collection with 3599 features and 18 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 0.782811 ymin: 6.240975 xmax: 3.821388 ymax: 12.40508
Geodetic CRS: WGS 84
First 10 features:
ID gid_0 gid_1 gid_2 gid_3 name_1 name_2 name_3
1 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
2 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
3 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
4 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
5 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
6 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
7 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
8 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
9 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
10 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
type_3 engtype_3 id_gns full_name full_nm_nd sort_name desig_cd
1 Arrondissement Borough 2058 Gano Gano GANO PPL
2 Arrondissement Borough 1426 Kokabo Kokabo KOKABO PPL
3 Arrondissement Borough 973 Wagui Wagui WAGUI PPL
4 Arrondissement Borough 3392 Atabénou Atabenou ATABENOU PPL
5 Arrondissement Borough 2726 Komon Komon KOMON PPL
6 Arrondissement Borough 2263 Dérou Derou DEROU PPL
7 Arrondissement Borough 3386 Banikoara Banikoara BANIKOARA PPL
8 Arrondissement Borough 453 Toké-Banta Toke-Banta TOKEBANTA PPL
9 Arrondissement Borough 1873 Gomparou Gomparou GOMPAROU PPL
10 Arrondissement Borough 581 Sonworé Sonwore SONWORE PPL
adm1 lat_y lon_x geometry
1 BJ-000 11.24685 2.414291 POINT (2.414291 11.24685)
2 BJ-000 11.28007 2.427724 POINT (2.427724 11.28007)
3 BJ-000 11.28146 2.400951 POINT (2.400951 11.28146)
4 BJ-000 11.28761 2.381083 POINT (2.381083 11.28761)
5 BJ-000 11.28951 2.423543 POINT (2.423543 11.28951)
6 BJ-000 11.29832 2.401679 POINT (2.401679 11.29832)
7 BJ-AL 11.29845 2.438561 POINT (2.438561 11.29845)
8 BJ-000 11.30770 2.414829 POINT (2.414829 11.3077)
9 BJ-AL 11.33134 2.441141 POINT (2.441141 11.33134)
10 BJ-000 11.35827 2.422149 POINT (2.422149 11.35827)
10.1 Calculating distance
ave_dist_rivers_points <- ave_dist_2shp(country_ext=gadm3,target=points_gns,shape_input=rivers_fao,UTM=31,res_raster=10000,name_dist="rivers")[1] "UTM is: 31"
ave_dist_rivers_points[1:10, ] #print first 10 elements ID gid_0 gid_1 gid_2 gid_3 name_1 name_2 name_3
1 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
2 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
3 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
4 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
5 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
6 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
7 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
8 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
9 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
10 1 BEN BEN.1_1 BEN.1.1_1 BEN.1.1.1_1 Alibori Banikoara Banikoara
type_3 engtype_3 id_gns full_name full_nm_nd sort_name desig_cd
1 Arrondissement Borough 2058 Gano Gano GANO PPL
2 Arrondissement Borough 1426 Kokabo Kokabo KOKABO PPL
3 Arrondissement Borough 973 Wagui Wagui WAGUI PPL
4 Arrondissement Borough 3392 Atabénou Atabenou ATABENOU PPL
5 Arrondissement Borough 2726 Komon Komon KOMON PPL
6 Arrondissement Borough 2263 Dérou Derou DEROU PPL
7 Arrondissement Borough 3386 Banikoara Banikoara BANIKOARA PPL
8 Arrondissement Borough 453 Toké-Banta Toke-Banta TOKEBANTA PPL
9 Arrondissement Borough 1873 Gomparou Gomparou GOMPAROU PPL
10 Arrondissement Borough 581 Sonworé Sonwore SONWORE PPL
adm1 lat_y lon_x ave_dist_rivers
1 BJ-000 11.24685 2.414291 5455.917
2 BJ-000 11.28007 2.427724 5483.610
3 BJ-000 11.28146 2.400951 4977.966
4 BJ-000 11.28761 2.381083 5037.752
5 BJ-000 11.28951 2.423543 5582.816
6 BJ-000 11.29832 2.401679 5329.104
7 BJ-AL 11.29845 2.438561 6046.705
8 BJ-000 11.30770 2.414829 5767.787
9 BJ-AL 11.33134 2.441141 6715.481
10 BJ-000 11.35827 2.422149 6887.261
10.1.1 Comparing results by polygons and by points
p1 <- hist(ave_dist_rivers_points$ave_dist_rivers)p2 <- hist(ave_dist_rivers$ave_dist_rivers)plot(p1, col=rgb(0,0,1,1/4)) # first histogram
plot(p2, col=rgb(1,0,0,1/4), add=T) # add
legend("topright", legend=c( "ave_dist_rivers_points","ave_dist_rivers"),
fill=c(rgb(0,0,1,1/4), rgb(1,0,0,1/4)), border=NA)